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ABSTRACT 



Wc investigate the column density distribution function of neutral hydro- 
gen at rcdshift z = 3 using a cosmological simulation of galaxy formation from 
the Overwhelmingly Large Simulations (OWLS) project. The base simulation 
includes gravity, hydrodynamics, star formation, supernovae feedback, stellar 
winds, chemo dynamics, and element-by-element cooling in the presence of a uni- 
form UV background. Self-shielding and formation of molecular hydrogen are 
treated in post-processing, without introducing any free parameters, using an ac- 
curate reverse ray-tracing algorithm and an empirical relation between gas pres- 
sure and molecular mass fraction. The simulation reproduces the observed z — 3 
abundance of Ly-a forest, Lyman Limit, and Damped Lj-a Hi absorption sys- 
tems probed by quasar sight lines over ten orders of magnitude in column density. 
Self-shielding flattens the column density distribution for A^hi > 10^^ cm~^, while 
the transition to fully neutral gas and conversion of Hi to H2 steepen it around 
column densities of A^hi = 10^°-^cm~^ and A^hi = 10^^-^cm~^, respectively. 

Subject headings: Methods: numerical — Quasars: absorption lines — Galaxies: 
formation — intergalactic medium — large-scale structure of Universe 



1. Introduction 



Ground-based spectroscopic observations targeting quasars are excellent probes of 2; > 1.7 
neutral hydrogen (e.g. |Rauch|[T998| [Wolfe et"aL][2005| . The Sloan Digital Sky Survey (SDSS) 
has produced approximately 1.5 x 10"^ moderate resolution quasar spectra (Abazajian et al. 



2009). These spectra provide ample data on Hi absorption lines with column densities 
A^Hi > 10^°-^ cm~^, so called Damped Ly-a systems (DLAs) (jProchaska k Wolfe]|2009[ |No- 



terdaeme et al. 2009). Lines with A^hi < 10 cm , the so called Ly-a forest, are best 



discovered in high-resolution spectra of bright quasars (e.g. [Kim et aL 2002). Lines with 
intermediate column densities, Lyman Limit Systems (LLSs), lie on the fiat part of the 
curve of growth, which complicates the determination of their column densities. Traditional 
methods of measuring A^hi in DLAs can be applied to high-resolution spectra for lines with 
Ahi > 10^^ cm~^ when damping wings begin to appear (e.g. Peroux et al.||2005l lO'Meara 
et al. 2007). Progress on the most difficult lines with 10^^'^ cm~^ 



< A^HT < 10^^ cm ^ has re- 



cently been made by Prochaska et al. (2010) by combining independent measurements of the 



Lyman limit mean free path and integral constraints over the column density distribution. 

Combining the observations above, one can determine the Hi column density distribu- 
tion function /(A^hi;-^), i-e. the number of lines per unit column density dN^, per unit 
absorption distance dX, at redshifts ^ ~ 3 from A^hi = 10^^ cm~^ to A^hi = 10^^ cm~^. Early 
determinations of /(A^hij-z) at these redshifts were reasonably well described by a single 
power law, /(A'hi, z) oc A'hi', with r] = 1.5 (Tytler 1987). As the quality of observations im- 



proved, this was no longer the case. Petitjean et al. (1993), showed that a single power law 
and a double power law with a break at A^hi = 10^^ cm~^ both failed Kolmogorov-Smirnov 



tests at the 99% confidence level. The most recent observations are well fit by a series of six 
power laws which intersect at A^m = {10^^•^ 10^^-3, 10^^ °, lO^o-^, lO^^-^S} cm'^ ([Prochaska 



et al. 2010). 



Attempts to explain the shape and normalization of /(Ahij -2) in a cosmological context 
have typically focused on sub sets of the full column density range. Analytic (e.g. Schaye 



2001a), semi-analytic (e.g. Bi & Davidsen 1997), and numerical (e.g. Theuns et al. 1998b||a ) 



models were instrumental in identifying the Ly-a forest lines with the diffuse, photo-ionized, 
intergalactic medium. Numerical work has also played a large role in determining properties 



of higher column density systems (e.g. 


Katz et al.||1996; Gardner et al.||1997; Haehnelt et al. 


1998 Cen et al.|2003| Nagamine et al. 


2004 Razoumov et al. 2006 


K^ohler & Gnedin 


2007 


Pontzen et al. 


2008 


Tescari et al. 2009 Hong et al. 2010 Cen 2010 


Nagamine et al. 


2010 


McQuinn et al 


2011 


)• 



Although self-shielding is crucial for modelling optically thick absorbers, only Razoumov 



et al. (2006), Kohler & Gnedin (2007), Pontzen et al. (2008), and McQuinn et al. (2011) have 
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used 3-D radiative transfer to calculate the attenuation of the UV background. Additionally, 



conversion of Hi to H2 is thought to determine the high end cut off in /(A^hi,^) (Schaye 



2001b Krumholz et al. 2009), yet only Cen (2010) included this process when modelling Hi 



absorption. We present a cosmological simulation of structure formation, to which we have 
applied a radiative transfer self-shielding calculation and a prescription for the conversion of 
Hi to H2 without introducing any free parameters. We show that this simulation reproduces 
observational determinations of /(iVni, z) around z = 3 over the entire range in column 
density. In addition, we determine the typical neutral fractions and total hydrogen number 
densities for Hi absorbers as a function of column density A^hi- 



2. Methodology 



We focus on model REF_WMAP7_L025N512 from the Overwhelmingly Large Simula- 
tions (OWLS) project ( |Schaye et ar][20To| ), which is identical to REF_L025N512 except that 
it was run using WMAP7 cosmological parameters. This simulation was performed with a 



modified version of the Smoothed Particle Hydrodynamics (SPH) code GADGET (Springel 



2005), and includes "sub-grid" models for star formation (Schaye & Dalla Vecchia 2008) 



chemodynamics (Wiersma et al. 2009a), galactic winds (Dalla Vecchia & Schaye 2008), and 



element-by-element cooling in the presence of a uniform UV background (Wiersma et al. 



2009b). Gas in the interstellar medium (ISM) at densities above n* = 0.1 cm~^ is assumed 
to be multi-phase and star-forming. This is modelled by imposing a polytropic equation of 
state (EoS) of the form P = P*(njj/n* )^/^. Because surface density and pressure are directly 
related in self-gravitating systems, the Kennicutt-Schmidt star formation law can be rewrit- 



ten as a pressure law (Schaye & Dalla Vecchia 2008). The observed Kennicutt-Schmidt law 



can then be used to determine a star formation rate in each gas particle. 

The simulation contains 2 x 512'^ particles in a periodic cube of size 25 comoving 
The cosmological parameters used are, {flra = 0.272, fib = 0.0455, Qa = 0.728, Ug = 
0.81, Us = 0.967, h = 0.704} (Komatsu et al. 2011). The mass resolution is mb = 1.47 x 
lO^/i~^M0 and mdm = 7.32 x 10^/i~^Mq for baryonic and dark matter particles, respectively. 
The equivalent Plummer gravitational softening length is e{z) = 1.95/(1 + z) proper h~^kpc 
at high z but is not allowed to exceed 0.5 proper h~\pc, a value reached at z = 2.91. 



REF.WMAP7.L025N512 included the Haardt fc Madau (2001) (HMOl) optically thin 



UV background from quasars and galaxies, but we apply a self-shielding correction in post- 
processing as follows. For each particle, we trace rays out to a distance Zray along Aray 



directions defined using the HEALPix algorithm (Gorski et al. 2005). We then compute 



frequency dependent optical depths along each ray and integrate over the HMOl spectrum 
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Fig. 1. — Hi column density distribution function, f{Niii,z), at z ~ 3; simulation results 
are shown as curves and observational data as symbols. The low Nm curve is obtained using 
mock spectra fitted with VPFIT. Self-shielding and H2 are unimportant in this range. The 
high iVni curve is obtained by projecting the simulation box onto a plane and includes self- 
shielding and H2. The gap around A^hi ~ 10^^ cm~^ separates low and high A^hi- Poisson 
errors on the simulation curves are always smaller than their thickness. We also show high- 
resolution observations of the Ly-a forest (Kim et al. 2002, "Kim02"), LLSs (Peroux et al. 
2005, "Per05"; O'Meara et al. 2007, "Ome07"), analysis of SDSS DLA data (Noterdaeme 
et al. 2009, "NPLS09"), and power law constraints (Prochaska et al. 2010, "POWIO", open 
circles are spaced arbitrarily along power law segments and do not represent A^hi bins or 
errors). 
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to calculate a self-shielded photo-ionization rate, r^*^''^, as opposed to the optically thin 
rate, T^^™ = rthin^]^g-i2 g-i _ I iQ This characterizes each particle with an effective 
optical depth Tes = — ln(r^'^^'^/r*^"^). We then use T^^^'^ to calculate a new neutral fraction, 
Xjjj = rijjj/njj, for each particle using an analytic equilibrium solution. We continue to 
loop over the particles until the neutral fractions converge. We obtain converged results for 
/(iVni, z) using Zj-ay = 100 proper kpc and A^j-ay = 12. Our self-shielding algorithm will be 
discussed in detail elsewhere (Altay et al., in prep.). 

In the OWLS snapshots, the temperature stored for gas particles on the polytropic 
star forming equation of state is simply a measure of the imposed effective pressure. When 
calculating coUisional ionization and recombination rates, we set the temperature of these 
particles to Tjsm = 10^ K. This temperature is typical of the warm-neutral medium phase of 
the ISM but our results do not change if we use lower values. We use case A (B) recombina- 
tion rates for particles with Tcs < (>)1. In addition, the optically thin approximation used 
in the hydrodynamic simulation leads to artificial photo-heating by the UV background in 
self-shielded particles. To compensate for this, we enforce a temperature ceiling of Tshid = 10^ 
K in those particles that become self-shielded {i.e. attain t^s > 1). In Figure 2 we show the 
effect of this temperature correction. 

For conversion of atomic hydrogen to molecules, we adopt a prescription based on ob- 



servations by Blitz & Rosolowsky (2006) of 14 local spiral galaxies to form an H2 fraction- 
pressure relation. Their sample includes various morphological types and spans a fac- 
tor of five in mean metallicity. They obtain a power law scaling of the molecular frac- 
tion, -Rmoi = SHa/^Hi, with the galactic mid-plane pressure, Rmoi = (-Pext/-Po)°, with 
a = 0.92 and Po/^b = 3.5 x 10^ cm~^ K. Applying this relation to the simulated ISM 
yields /h, = [l + A {njniy"]-' with A = (P./Pq)"", and /3 = «7eff. 

The Hi column density distribution function, 

„,,^ , (fn (Pn dz 

f{N-iii,z) = 



dNuidX dNmdzdX ' ' 

is defined as the number of absorption lines n, per unit column density dNfn, per unit 
absorption distance dX. The latter is related to redshift path dz as dX/dz = Ho{l + 
zy/H{z), where H{z) is the Hubble parameter ( Bahcall fc Peebles|l969 ). In the comparisons 



below, we scale f{Ni{i,z) reported by various observers to the cosmology assumed in our 
simulation. 

The simulated f{Nm,z) below A^hi = lO^^cm"^ is computed by generating 1000 mock 
spectra through each snapshot. We then apply instrumental broadening with FWHM 6.6 km 
s~^, add Gaussian noise such that we have a signal-to- noise ratio of 50 in the continuum, and 
fit the mock spectra using VPFIT ( Carswell et al.|[l987 ); see Theuns et al. (1998b) for more 
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details. To obtain /(A^m? -2) for the rarer systems with N^ii > lO^^cm"^, we project all 512^ 
gas particles along the z-axis onto a grid with 16,384^ pixels using gaussian approximations 
to their SPH smoothing kernels. This leads to hypothetical lines of sight with a transverse 
spacing of 381 proper h~^pc or about 3/4 the gravitational softening length at 2; = 3. We 
have verified that our results are converged with respect to the projected grid resolution. 
For systems with A^hi > 10^^'^ cm~^ and redshifts z < 4.4, the rate of incidence per unit 



absorption distance, 1{X), is observed to be less than one (Prochaska et al. 2010). The 
absorption distance for a single sight-line through our box at z = 3 is AXi = 0.133 and so 
we expect, on average, much less than one system per sight-line. Therefore, the contribution 
to the total column density in projected pixels with A^hi > lO^'''^ cm~^, for the vast majority 
of cases, is dominated by the single absorption system in the line of sight. Curves in Figure [T] 
are labelled either "VPFIT" or "Projected" , depending on the method used, all others were 
calculated using projections. 

Table 1 hsts the iVni bins, absorption lines per bin, and total absorption distance used 
for the low and high A^hi analyses of our fiducial model. The (25/i~^Mpc)^ volume searched 
for absorbers contains ~ 39, 000 friends of friends dark matter halos with masses above 
7.32 X lO^/i~^M0 and yields ^ 2 x 10^ lines of sight containing DLAs. The size of this 
data set obviates the need to re-weight a limited sample of absorbers using an analytic mass 



function as in Gardner et al. (1997) or Pontzen et al. (2008). 



Results 



3.1. Full Range 



1022 cm-2. 



In Figure[T| our fiducial model /(A^'hi, z) is plotted at ^ = 3 from A^hi = 10^^- 
The analysis using VPFIT in the Ly-a forest range, where self-shielding and H2 are not 
important, joins smoothly onto the projection analysis at A^hi > 10^^ cm"^. The model is 



compared to high-resolution observations of the Ly-a forest (Kim et al. 2002) and LLSs 



dPeroux et aL]|2005[ [O'Meara et 'aL]|2007D , DLA statistics from the SDSS ( |Noterdaeme et al. 



(Peroux et al. 


2005 


2009 


), and a series 



Both our model /(A^hi, z) and the observations display a characteristic flattening above 
the transition to Lyman Limit Systems at A^hi = 10^^'^ cm"^^ and a steepening beginning 
around the DLA transition, A^hi > lO^^-^ cm"^. To quantify this model's goodness of fit to 
the data, we calculate per degree of freedom between the model and the three largest data 
sets using the error bars reported by the observers. The corresponding poisson error bars for 
our model are smaller than the thickness of the curves shown in Figure [Tj The results are 
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Table 1: Simulation line list. 



VPFIT Projection 



1000 X AXi 


= 133.1 






16, 


384^ X AXi 


= 3.574 X 10^ 




A log 




^ of lines 


A log 
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1.26 for the Ly-a forest data from Kim et al. (2002) and 1.26 and 2.24 for DLA data from 



Prochaska & Wolfe (2009) and Noterdaeme et al. (2009), respectively. Lower normalizations 



of the UV Background as found in Haardt & Madau (2011) and shown in Figure 2, would 
improve these fits. 



3.2. LLSs and DLAs 

In the left panel of Figure [2] we show models in which the amplitude of the UV back- 
ground was varied by factors of 3, a model with no self-shielding, and our fiducial model. 
Although systems with A^hi > lO^'''^ cm~^ are optically thick to photons at the Lyman limit, 
models with and without self-shielding (at the fiducial UV background normalization) don't 
diverge until A'hi = lO^^cm"^. This is because higher energy photons with smaller cross- 
sections for absorption penetrate the clouds. Between A^hi = 10^'' cm~^ and A^hi = 10^^ cm~^, 
the model with self-shielding predicts fewer lines, because systems are moved to higher col- 
umn densities in the self-shielded model. 

Above Ahi = 10^^ cm^^, the model that neglects self-shielding stays on the Ly-a forest 
power law until it steepens around Ahi = 10^^'^ cm~^ due to the formation of molecules. 
The other models flatten due to self-shielding and then steepen due to both the formation 
of molecules and the saturation of the neutral fraction. The flattening of f{Nm,z) is a 
hallmark of self-shielding and was also found in the original numerical work of |Katz et al. 



(1996) and in the analytic work of Zheng & Miralda-Escude (2002). Changes in the UV 



background normalization by factors of three result in constant shifts of f{Ni{i,z) until the 
gas is completely shielded around A'hi = 10^^'^ cm~^. This normalization adjustment is larger 



than any of the uncertainties claimed in recent work (e.g. Faucher-Giguere et al. 2008) 



3.3. DLAs 



In the right panel of Figure |2j we isolate the effects of H2, and the photo- heating of 
self-shielded gas. The models with H2 approach a vertical asymptote just above A'hi = 
IQ22.0 t^IiHq ^]2e model without H2 predicts the existence of systems out to A^hi = 

]^g24.5 (3)^-2 although at such low abundance, that less than one would have been discovered 
in the SDSS. 

The introduction of H2 produces a steepening of /(A'hi,^;) around A'hi = lO^^'^cm"^. 
Such a transition, suggested theoretically in Schaye| (2001b), has been observed at 2; = 
using CO maps as a tracer for H2 (e.g. Zwaan & Prochaska 2006). This feature coincides 
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Fig. 2. — f{Ni{i, z) - LLS and DLA range. In the left panel, we vary the amplitude of 
the UV background and show the impact of neglecting self-shielding. In the right panel, we 
isolate the effects of H2 and show a model in which we have not lowered the temperature in 
self-shielded particles (w/o Tghid)- On top of each panel, we show the ratio of each model to 
our fiducial model (solid red curve), which includes self-shielding and II2. The observational 
data are a subset of those in Figure [T] plus SDSS analysis from Prochaska & Wolfe 2009, 
"PW09", in the right panel. Self-shielding becomes important for A^hi > 10^^ cm~^ leading 
to a flattening of f{Niii,z). Cooling the self-shielded gas yields a constant offset while H2 
becomes important above column densities of iVni > lO^^'^. 



- 10 - 



with the break in the double power law commonly used to fit /(A^m, z) in the DLA column 
density range (Prochaska & Wolfe 2009 Noterdaeme et al. 2009[ ), suggesting a relationship 
between the two. At DLA column densities, ionizing radiation from local sources may play 



a role (Schaye 2006). We have not included these sources in our self-shielding model, but 



Nagamine et al. (2010) have recently shown that /(A'hi, z) changes by less than 0.1 dex when 
local sources are included. 

Because the UV background suppresses cooling, the temperature recorded in the OWLS 
snapshots for particles that are identified as self-shielded in post processing is an over- 
estimate. To compensate for this, we enforce a temperature ceiling of Tshid = lO^K in 
self-shielded particles in our fiducial model. The curve labelled "w/o Tshid" shows a model 
in which we have not performed this temperature adjustment. Because the temperature 
dependence of the collisional equilibrium neutral fraction is very small below lO'^K, the two 
temperature models should bracket the neutral fractions one would expect from a more 
accurate treatment of the temperature. The difference between these two models is about a 
factor of ten smaller than the difference between the optically thin and self-shielded models 
but can be on the order of the observational one sigma error bars around the DLA threshold 
where the data are most abundant. We plan to explore hydrodynamic simulations that 
include self-shielding in future work. 



3.4. Physical properties of high A^hi absorbers 

Neutral hydrogen mass weighted values for the neutral fraction, x^^ = n^Jn^, and total 
hydrogen number density, = n^^ + n^^^ + '2n^^ are plotted as a function of A^hi in Figure 
[3| The effects of self-shielding produce a steep deviation from the optically thin power law in 
Xjjj above A^hi = lO^^cm"^. As the UV Background normalization is reduced, and as higher 
temperatures are used, the deviation becomes smaller. The median x^^ at A^hi = 10^^ cm~^ 
in our fiducial model is 0.3, however there is a large spread in the data in this column density 
range. It begins to drop around A^hi = 10^^ cm~^ due to the formation of H2. Systems above 
A^Hi = 10^^ cm~^ have lost much of their atomic Hydrogen to molecules, however the H2 
likely has a small covering fraction. 

The median flattens around the beginning of the LLS range, A'hi = 10^'''^cm~^, 
to approximately 2 x 10^^ cm^^ where it remains roughly constant until the start of the 
DLA range, A'hi = 10^°'^cm~^. Above this column density, the gas is fully neutral (see left 
panel) causing to rise steeply with A^hi and /(A^hi,^) to steepen (see Figure [2]). Above 
A^Hi = 10^^ cm~^, the medians for models which include H2 are steeper than linear due to 
the formation of molecules. The normalization of the UV Background and the treatment of 
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Fig. 3. — Left panel: n^j weighted neutral fraction, Xjjj(iVHi). The red sohd hne indicates 
median values in A^hi bins for the fiducial model which includes self-shielding and H2. The 
contours represent 68%, 95%, and 99% of the data about this median in each bin. Also shown 
are median values for the models shown in Figure 2 and a model with lower UV background 
normalization and no temperature adjustment for self-shielded gas. Right panel: As previous, 
but for the rij^j weighted total hydrogen number density 



^Hi + ^Hii + 2^H2 • show 



the predictions of the analytic, optically thin model of Schaye (2001a). H2 begins to reduce 

1020.3 , 



Xjjj around the DLA threshold 



HI 



cm 

compared to the optically thin case between 10^*^cm~^ < A'hi < 10 



and self-shielding flattens the median 



20.5 2_ 



temperature can change the LLS characteristic density by half a decade. For the optically 
thin case we find excellent agreement with the corresponding prediction in Schaye (2001a). 



4. Conclusions 

We have used a hydrodynamic simulation of galaxy formation together with an accurate 
ray-tracing treatment of self-shielding from the UV background and an empirical prescription 
for H2 formation, to compute the z ^ 3 Hi column density distribution function. We find 
agreement between the reference OWLS model and the entire column density range probed 
by observations (lO^^cm"^ < A'^hi < lO^^cm"^). We have shown that f{Nm,z) fiattens 
above A^'hi = 10^^ cm~^ due to self-shielding, and steepens around A'hi = 10^°'^ cm~^ and 
A^Hi = 10^^'^ cm~^ due to the absorbing gas becoming fully neutral, and the transition from 
atomic to molecular hydrogen, respectively. In future work, we will examine the systems 
causing this absorption in greater detail and repeat these analyses on a large sample of 
OWLS models. 
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